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The localization of vibrational energy, induced by the modulational instability of the 
Brillouin-zone-boundary mode in a chain of classical anharmonic oscillators with finite ini- 
tial energy density, is studied within a continuum theory. We describe the initial localization 
stage as a gas of envelope solitons and explain their merging, eventually leading to a single 
localized object containing a macroscopic fraction of the total energy of the lattice. The initial- 
energy-density dependences of all characteristic time scales of the soliton formation and merging 
are described analytically. Spatial power spectra are computed and used for the quantitative 
explanation of the numerical results. 
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I. INTRODUCTION 



Despite of a large number of studies on localization of energy in pure anharmonic lattices 0-H] , the mechanism 
of their thermal generation is still under debate. This is clearly a central issue being directly related with the 
possibility of observing signature of such a class of excitations in experiments on real crystals. The problem 
of nanoscale energy localization has turned out to be rather important also in view of its possible application 
for different physical and biological systems. The theoretical studies at finite energy per particle dealt mainly 
with one-dimensional models belonging both to the Fermi-Pasta- Ulam (FPU) class of model [g]]6) (with acoustic 
spectrum in the harmonic limit) and to the Klein-Gordon one M (with optical-phonon-type spectrum in the 
• i-H , harmonic limit). Some attention has been also recently payed to the discretized version of the nonlinear 
^^ ' Schrodinger equation in one dimension |jjj£)| . 

Of special interest is the numerically observed fact that excitation of short-wavelength extended modes of 
the chain and their instability with respect to perturbations with nonvanishing wavenumber leads generically to 
creation of localized states. The phenomenology emerging from the numerical studies |g|| can be summarized 
in three main stages. In the first one, after the development of the modulational instability, energy localizes 
along the chain in the form of an array of almost standing "bumps" . After a given time (that increase upon 
decreasing the temperature) those begin to merge diminishing their number and localization lengths. This 
second stage (which is reminiscent of the phenomenon of "soliton turbulence" ||10| ) ends with the emergence of 
a single localized and moving object (a "breather"), carrying along a macroscopic fraction of the total energy. 
In the third stage, the breather continues to slow down radiating its energy towards the background until it 
eventually disappears and the equipartition (thermal equilibrium) state is attained. 

The aim of this paper is to give a consistent quantitative explanation of such behaviour, in particular of 
the apparent paradoxical tendence of the system to increase localization as time goes by. Namely, we would 
like to clarify the role of the finite vibrational energy density on the formation of strongly localized envelope 
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solitons and breathers and their temporal evolution in the anharmonic FPU chain. This is accomplished 
by means of a suitable continuum approximation for the lattice dynamics. The approach in terms of the 
envelope-function equation including higher-order nonlinear terms has been already introduced to describe 
short-wavelength excitations in the FPU chain |ll[ and is reviewed in Section II. It is a useful framework both 
for computing localized solution in the form of envelope soliton (ES) and breather (section III), and for describing 
the modulational instability (Section IV). Furthermore, going beyond the usual integrable approximation |12(| 
is of crucial importance to explain the localization process (Sections V and VI). We described analytically the 
initial-energy-density dependence of all characteristic time scales of the envelope soliton formation and merging. 
Particular emphasis has been also given to the analytic computation and numerical measurements of the spatial 
power spectra ("structure functions") that, at least in principle, are of importance for real experiments. Indeed, 
they turn out to be extremely useful for the physical interpretation of the simulation results (see, e.g., Sections 
V and VII). Finally, a brief discussion of the finite-temperature effects is reported in Section VIII. 

II. ENVELOPE-FUNCTION EQUATION FOR THE FPU MODEL 

We start from the classical FPU Hamiltonian of a monatomic chain of N particles of mass m with anharmonic 
nearest- neighbor potential of order r > 3: 
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where u n is the (real scalar) displacement of the n— th particle from its equilibrium position and K 1 are the 
harmonic (7 = 2) and anharmonic force constants. The lowest order potential describing the anharmonicity of 
the longitudinal or pure transverse motion in a centrosymmetric ID lattice corresponds, respectively, to r = 3 
or r = 4. From the Hamiltonian (uh we obtain the following equations of motion: 
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In the following, periodic boundary conditions are always assumed. 

In order to describe the short-wavelength excitations of the chain with wavenumber k ~ tt/cl (a being the 
lattice spacing), i.e. near the Brillouin-zone boundary, it is convenient to introduce an envelope- function of 
the displacements f n = (— l) n u n , /„ = f(x/a). The latter is assumed to be slowly varying on the interatomic 
scale: af x -C /• Substituting in Eq. (||) the expansion of the differences u n ±\ — u n up to the fourth order, one 
can obtain a nonlinear partial differential equation for / fllj]. For instance, in the case of the purely quartic 
anharmonicity (i.e. the so-called /3-FPU model with only Ki and K4 different from zero) such an equation reads 



mf + AK 2 f + K 2 a 2 f xx + -^K 2 a 4 f x 



+ 16K 4 f + 6K 4 a 2 f(f) xx = 



(3) 



where the subscript x stands for the spatial derivative of the corresponding order. Equation, similar to Eq. (H), 
can be also obtained for the cubic anharmonicity (a-FPU model with 7 = 2 and 7 = 3) and, more generally, 
for the FPU chain (p with the anharmonic potential of the arbitrary (even or odd) order r (u| . 

Due to the last term in l.h.s. of Eq. (g), the latter substantially differs from the conventional nonlinear Klein- 
Gordon and its nonrelativistic version - the nonlinear Schrodinger equation. Actually, Eq. (O) can be reduced to 
the latter only in the limit of relatively small oscillation amplitudes (see Section IV below) . Eq. (g) differs also 
from the well-known Korteweg-de Vries equation which describes the evolution of the long-wavelength (i.e. the 
"acoustic-phonon-type" ) excitations and can be derived from the lattice Hamiltonian (nl) in the limit au x <C 1 
(see, e.g., p)). 

In what follows, we will use Eq. (0) for the adequate analytic description of all the main stages of modulation- 
instability-induced dynamic evolution of the hard (with A" 4 > 0) quartic FPU chain. For simplicity of notation, 
we set K 2 , K4, a, and m to unity. With such a choice, the (adimensional) energy density (i.e. the energy per 
particle) e is the only relevant physical parameter if one does not account for the finite-temperature effects (see 
Section VIII below). Moreover, we will always deal with the (classical) low-temperature limit that in our units 
corresponds toe< 1. 



III. LOCALIZED SOLUTIONS: ENVELOPE SOLITONS AND BREATHERS 

As a first step, we determine the form of the single ES (or periodic array of solitons) by looking for solutions of 
Eq. (|3J) in the monochromatic form f(x, t) = 4>{x) cos tot. The frequency ui can be computed by first substituting 
such an Ansatz in (||) and invoking then the so-called Rotating Wave Approximation (RWA), namely to replace 
cos 3 cot by j cos cut. It is worth mentioning that such an approximation is expected to hold for the considered 
short- wavelength vibrations due to the weakness of nonresonant interaction between the mode with fundamental 
frequency and its third harmonic. By further neglecting the fourth-order dispersive term in Eq. (|3|) we find its 
first integral as Q: 

(4 - lj 2 )<J) 2 + (1 + 90 2 )Oz) 2 + 60 4 = const . (4) 

The generic solution of Eq. (0) for a nonvanishing constant corresponds to a periodic array of ESs, and can 
be conveniently parametrized by the two extreme values of the amplitude 4>min and 4> ma x (so that 4> m i n < ^ < 
4>max)- Those are simply determined by the conditions that 4> x vanishes for <j) = 4> m in and 4> = 4>max- From 
those two conditions we derive immediately the frequency of the envelope-soliton array as a function of the 
extreme values of its amplitude: 

Lu 2 = 4 + 6{tf mm + tf max ) . (5) 

Then from Eqs. (§]|3|) we can determine the desired solution as 

t>Y_ A<f> 2 m a X -4> 2 )(<t> 2 -4> 2 min) , fi x 

dx) 1 + 90 2 ' { ' 

Eq. (f|) can be solved by quadratures yielding the following implicit form for 4>{x): 
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The single ES corresponds to const = in Eq. (Q) and to cf> m i„ = in Eq. (0). Although a general explicit 
form of Eq. (m) is not feasible, we can at least discuss the two limiting cases of small and large amplitudes. In 
the small-amplitude limit (9<f>ma,x *C 1), the solution is given by (we include explicitely the time dependence) 

t <„ +\ A COs(uj s t + lf>) 

fs{x,t) = 4> m ax — — 7=- , (8) 

cosh[(a; - x o )0 mQX v6J 

where the frequency as given by Eq. (||) with (f> m in — and is approximately cu s = 2 + 3(/> 2 na:E /2, whereas <p is 
the initial phase. 

In the large- amplitude limit (90 2 rlQ:! . ^> 1), it follows from Eq. (0) that the ES acquires a sinusoidal shape 
with a short amplitude-independent localization length |11| 

7T 

fb(x,t) = (f> max cos(uJbt + ip)cos[h(x — xq)] for fo\x - xq\ < — (9) 

and with /b sa outside of the above specified domain. The effective wavenumber is fcb = v/2/3 and lo^ 
given by Eq. (g) with <f> m in — 0. This is rather accurate approximation for the exact breather solution of 
the discrete FPU model (see Ref. [|3[ for a detailed comparison). Discrepancies have to be expected because, 
strictly speaking, the continuum approach is not completely justified in this short-wavelength limit. Indeed, the 
analysis of the corresponding discrete-lattice equations predicts |14| the existence of intrinsic localized modes 
(discrete breathers) with (approximately) sinusoidal envelope (g) but with k b = 7r/3 and lu 2 = 3 + (81/16)<f> 2 nax . 
Having at hand the explicit solution (|8|), we can compute the spatial Fourier transform of the displacement 
field u associated to the single ES obtaining fj 



1 Notice that, by definition, the zero wavenumber for / corresponds to the wavenumbers ±7r for u and therefore u(k) 
J u{x)e lkx dx = / f{x)e lkx cos irxdx = \ [f{k -ir) + f(k + it)] . 



u s (k,t) = -cos(u s t + <p) [F s (k-w, x Q ) + F s (k + ir,x )} , (10) 

where we have defined the form factor 

F s (k lXo ) = . (11) 

Accordingly, the power spectrum (i.e. the square modulus of u s (k, t)) is exponentially localized around k = ±7r 
and its maximal value is independent of the amplitude 4> max . 

For later purposes, we compute also the total energy E s of the ES. To this aim, one can use the expression 
for the hamiltonian density h by the variation of which the equation (pft can be derived Jllj ] 

h = \(f? + 2/ 2 + 4/ 4 - i(l + 12 f 2 )(f x ) 2 + i(/ x ,) 2 , (12) 

(remember that we are working in the dimcnsionless units). Neglecting again the higher-order dispersive term 
and averaging h over time (in the RWA) and by means of (^) with const = 0, we find that the energy density 
h B associated with the ES is 2</> 2 + 3</> 4 . Using this expression together with (||), we finally find 

/+oo /-+oo A 

h s dx= / [2cf> 2 + 3q> 4 ] dx = -=ct> max (l + <f> 2 max ). (13) 
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In the leading approximation, the energy scales linearly with the amplitude [ p0[ . 
The Fourier transform of the breather solution (Bh is 

Ub(k,t) = -cos(w&i + <p)[F 6 (fc-7r,a;o) + i'&(fc + 7r,a:o)] , (14) 
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where we have defined 

F b (k,x ) = ( f >max e tkx ° 1 ^P^cos\^-\ . (15) 

Notice that this function is finite for k = k^. In particular, for kb — n/3 the Fourier transform of the breather 
solution has the following characteristic form: 
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The power spectrum has basically a sinusoidal shape with two absolute maxima at k — ±7r whose height is 
proportional to the square of the breather amplitude. 

IV. MODULATIONAL INSTABILITY 

In this section, we describe the modulational instability of the Brillouin-zone-boundary mode corresponding 
to the uniform envelope f{x,t) = fo(t) = Acosut. To compute the frequency u) one can resort again to the 
RWA mentioned in the previous section. In the small- amplitude limit 3 A 2 <C 1 one obtains u> = 2 + 3A 2 , which 
amounts to a correction to the zone-boundary frequency ui max = 2\/K^Jm = 2. 

Let us consider a spatially inhomogeneous small perturbation Sf oc cosfc m x, where the wavenumber k m of 
the envelope-function modulation is assumed to be much smaller than the Brillouin-zone width 7r, an hypotesis 
that will be confirmed a posteriori. Then, from the linearization of Eq. (0) we find the equation for 5f to be 

5f + [4 - k 2 n + 24 A 2 (I + cos 2u>t)]6f = . (17) 



To obtain the growth rate s at the parametric resonance, described by Eq. (|17|), we take a solution of this 
equation as (cf., e.g., fl5|0 . 



Sf = [a cos tot + bsmujt\e cos k m x , (18) 

where the unknown amplitudes a and b are assumed to be small with respect to A. Therefore, from Eqs. 
and rtlq) we obtain the dispersion equation for the growth rate s: 
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The maximal growth rate s max and the corresponding wavenumber k* m of the envelope-function modulational 
instability are thus given in the small-amplitude limit 3A 2 <C 1 by 



k* m = AV12 , s max =3A\ b=-a, (20) 

which, consistently, describes a longwavelength modulation with k* n /2 -C 1 |Q. In the considered limit, such a 
value of s max coincides with the one computed from the stability analysis of the discrete- lattice equations ]Pfl . 
It is also worth to mention that there is a minimal (threshold) amplitude A c for the modulational instability of 
the Brillouin-zone-boundary mode in a finite FPU chain of N particles, and A c ex 1/N for N 3> 1, see, e.g. 0. 
We neglect such finite-size effects in our continuum approach which implies, in particular, that the amplitude 



A in Eq. (20) is well above the threshold one for the considered (finite) chain. 

By means of Eq. (||) we can also describe the interesting phenomenon of the further appearance of spatial 
harmonics of the "primary" perturbation i.e. the components with wavenumbers nk^ and n — 0,2,3,.... To 
describe for instance the growth of the zeroth and second spatial harmonics, we assume Sf to be of the following 



form (cf. Eqs. (|l§) and fl20|)): 

Sf = a [coscot — sinujt] e Smaxt cosk^x + [ccosuit + dsincot] e 2Smaxt cos2fc* l ir 

+ [ecosujt + gsmtut] e 2Sm ^ 1 = Sf Ki + Sf 2k * m + Sf , (21) 

where the amplitudes c, d, e, and g are assumed to be much smaller than both the condensate amplitude A 
and the amplitude a of the primary perturbation. Using Eq. (0), we obtain the following evolution equations 
for 5f2k* and Sfo similar to Eq. (nTl) with a driving term — 48/o(<5/fc* ) 2 on the r.h.s.: 



Shk L + [4 - 4fc r *„ 2 + 24A 2 (1 + cos2ujt)]Sf 2kL = 

-12a 2 A[2cosujt~smujt]e 2Sma:ct cos2k^ n x , (22) 

Sfo + [4 + 24A 2 (1 + cos2ujt)]Sf = -12a 2 A[2coscj^ - sinwt]e 2s — * . (23) 

^From Eqs. (pi|)-(|23|) we find the amplitudes c, d, e, and g: 

a 2 , a 2 3a 2 .„„, 

c= — 7, d = e = -, q= . (24) 
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Therefore the zeroth 5f and second Sf 2 k* spatial harmonics grow exponentially with a rate which is two times 
larger than that of the primary instability Sf k * ■ Nevertheless, as their initial amplitudes are of second order 
with respect to the amplitude a of the primary instability, one expects that they will come into play only at a 
later stage. 

More generally, we can show that the growth rate of the n-th spatial harmonic Sf n k* of the primary instability 
is equal to ns max while its initial amplitude is proportional to the n-th power of the amplitude a. Furthermore, it 
can be shown that small contributions to the growth of the zeroth harmonic Sfo are induced by the higher-order 
harmonics. These arise from the terms like (5/2fc* n ) 2 /o, (<5/fc^) 2 ^/2fc* n etc. to the r.h.s of Eq. (p3[). 

The above scenario is illustrated in Fig. |l|, where the outcomes of a numerical simulation of the quartic 
FPU model are reported. The equations of motion were integrated with a third-order symplectic algorithm |18|] 
(microcanonical simulation) starting from the initial conditions 



u„(0) = 0, u n (0) = (-!)" V2e. (25) 



Actually, a small gaussian white noise (with amplitude <§; V2e) has also been added to the initial velocities in 
order to speed up the instability. The left panels of Fig. [j] show the evolution of the energy density (cf. Eq. 

(1)), 
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at subsequent times. The spatial power spectrum of the velocities 

N 9 4 

shows the growth of the perturbation with the wavenumbers tt — k^ (vertical line). Obviously a symmetric 
peak at — 7r + k* m is also present as well as the further transfer of the energy towards the spatial harmonics 
with wavenumbers nk^, n = 0, 2, 3, ..., according to the above described mechanism. We also checked that the 
typical time for developing the instability scales as A~ 2 oc e _1 as prescribed by Eq. (f20[). 

V. THE FIRST LOCALIZATION STAGE: A SOLITON GAS 

As already mentioned, the modulational instability produces a strong localization of the vibrational energy 
along the chain (see Fig. H). We wish now to describe such a set of localized objects. More precisely we will 
show that it is basically a gas of ESs and study some of its properties. 

Let us consider a set of N s non-overlapping ESs of the form <M) centered at the points Xi along the line and 
with initial phases ipi . For the sake of simplicity, let us also assume them to have approximatively the same 
amplitude. According to (Ru), the spatial Fourier transform of the corresponding displacement field will be 

1 Ns 
u(k,t) = -22cos(u} s t + ipi)[F a (k-Tr,xi) + F s (k + 7r,Xi)] . (28) 

i=i 



The velocity power spectrum Sk as defined by Eq. ( p,7\ j in the discrete case, can thus be approximated as (we 
use the explicit form of F s ) 

s k = ^(1 «(M) I 2 ) = ^f i cosh2[ !,(.-») I £e*«+«*->" I 2 
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where the angle brackets denote the averaging over a time interval larger than 2tt/uj s (that is performed to 
remove the oscillations as well as statistical fluctuations). In the last expression we neglected exponentially 
small overlap between the components of (|10|) localized respectively around k — tt and k = —tt. 

Finding a more explicit expression for 5^ requires of course some hypotesis on the statistical properties of 
the positions xi and phases ipi. For instance, in the limit of completely uncorrelated ESs one can replace the 
square modulus of each sum in Eq. (|2^) by N s . To take into account weak partial correlations, we introduce 
a phenomenological parameter p (p ~ 1) and replace the mentioned square modulus by pN s . Obviously, p = 1 
corresponds to a completely random case while the cases p > 1 or p < 1 describe, respectively, the arrays of the 
partially correlated or anticorrelated ESs. With such a definition 

_ tt 2 uj 2 s r 1 1 

Sk ~ PUs ^ { cosh 2 \ _(__! 1 + cosh 2r ^+«) I 1 ' m 
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where n s = N s /N is a density of the envelope solitons. 

To compare the latter expression with the numerically measured Sk, one has to relate u) s , n s and (\> ma x in 
Eq. (pfl) with the energy density e. As we deal with the limit e -C 1, one has that lo s ss 2 since A ss \/2e/2 (cf. 
Eq.(Ea) and the definition of the amplitude A). Moreover, the density n s will be given by the inverse wavelength 
of the modulational instability, namely 



This amounts to say that a single localized object emerges from each wavelength. A similar reasoning allows us 



to relate the total energy E s (which is basically proportional to (j> m ax, see again Eq. dl3|)) to e as 



E '=*i~***i\- (32) 

The phenomcnological parameter q, < q < 1, gives the fraction of the initial vibrational energy which is 
absorbed by the set of ESs (and actually accounts for the partial overlap between them). In particular, for 
e — > we expect the soliton gas to be more and more rarified and to contain almost all of the initial energy in 
the lattice so that q — > 1. 

A consequence of Eqs. (|3l]) and ( p2\ ) is that the ratio between the total energy of the ESs and their density 
is (almost) independent on e: 

^= q —- (33) 

Before proceeding further, we numerically checked that the above picture is correct by evaluating the number 
and energy of solitons as a function of the initial energy density. The results reported in Fig. ra show a very 



good agreement with both Eqs. (31) and (B3|). Notice that, as expected, the data are very well fitted with q = 1 
at least for e < 0.01. 

Finally, using Eqs. (|3(]), (|3l|), (p2|), and (|l3|), we obtain the desired expression for the velocity power spectrum 
as a function of e and the (unknown) parameters p and q only: 

Sk=P T6^ C oJ[^] + cosh 2 W] } ' (34) 

L gv6e L Qv6e J 

where lu s = 2. 

We obtained a very good fit of our numerical data, for energy density e = 0.0033, using Eq. (|34| ) with p = 0.50 
and q = 1.0 (see the right panels of Fig. ||). We interpret this by saying that the system is assimilated to a 
diluted gas of nonoverlapping and partially anticorrelated (i.e. out of phase) ESs. For larger initial energy 
density e = 0.033 the data are fitted by p = 1.12 and q — 0.78 (lower panels of Fig. g). This latter case 
correspond, as expected, to a more dense gas of partially overlapping ESs, which is characterized, in comparison 
with the previous one, mainly by pair correlations between neighbours. 

VI. THE SECOND LOCALIZATION STAGE 

To understand why a set of almost nonoverlapping small-amplitude ES finally merge into a large-amplitude 
one, it is first of all convenient to recast the second order (in time) equation (0) for the real envelope-function 
in two first order equations. For small-amplitude oscillations with a frequency close to the band-edge uj max = 2 
this is accomplished introducing the complex function ^(x^t) such that 

f(x,t) = i^(x s t)e-*-»"« + rfotjc*"""*] . (35) 

Substituting this expression into Eq. (S) and assuming that ip ^C Lu ma x4'} after some algebra one can obtain (in 
the RWA) the following coupled nonlinear-Schrodinger-type equations for -0 and ip* 

-iu> max ip + -ip xx + —Ipxxxx + 6 | ip | 2 ip 

+^(r^xx+ 1 i> x i 2 ) + i> 2 r xx + W%] = , 

iuJmaxlp* + ^xx + l^Plxxx + 6 | | 2 0* 

+|[2v>*(^: x + 1 i>x i 2 ) + r 2 ^xx + Wx 2 } = o . (36) 



Equations ( pq ) can be expressed in canonical form as 

i_sn _. i*_ 5 2L 

dtp dip 

where the Hamiltonian H = J Hdx has the following density H: 

ff = 4[|^l 2 -^IV^I 2 -6M 4 



+6[|^| 2 |^| 2 +i(^V: 2 + ^* 2 )]]. (37) 

One can ascertain that equations (B6J) admit the integrals of motion 

H=jHdx, V = - % - [[ipil)*-tp*tl> x ]dx , Af= J \ip\ 2 dx , (38) 

which are the energy, momentum and number of quanta, respectively. 

The inclusion of higher-order and nonlinear dispersive terms make equations (B6h nonintegrable (i.e. they 
do not have an infinite number of constants of motion). This simple fact, together with the conservation laws 
( p8| ) has deep consequences on the dynamics: it is in fact known |]1| that in this case the total number of 
solitons is not conserved. Actually, they can merge together enhancing their amplitudes and decreasing their 
number due to an exchange of both energy and number of quanta. Indeed, the Hamiltonian (|37|) describes an 
effective repulsive interaction between quasiparticles with negative (and amplitude-dependent) effective mass 
and it is therefore equivalent to a system of quasiparticles with positive effective mass and effective attractive 
interaction. As the solitons merge, the released binding energy is carried away by free waves and the "phononic" 
part of the spectrum is gradually populated. 

More precisely, the kinematics of the soliton-soliton interaction process requires an increase of the number of 
quanta in the stronger soliton and a decrease in the weaker one accompanied by the release of some binding 
energy by means of free plane waves that give also the recoil momentum to the solitons |L9| . Accordingly, 
the reverse process should include a triple collision (two solitons and a free wave) and therefore is much less 
probable than the direct one. It is therefore clear that eventually only one large-amplitude soliton survives on 
the background of small-amplitude free waves with a continuous frequency spectrum. This argument explains 
qualitatively why the purely dynamical process of soliton merging is time-irreversible. 

The above reasoning can be further carried on to give a more quantitative description of the phenomenology. 
For instance, it is possible to explain the observed fact that the characteristic time for complete soliton merging 
in the quartic FPU chain scales with the energy density as e~ 2 || (see also ||). To this aim, let us estimate 
the ES lifetime (or relaxation time) t ss due to the scattering mentioned above. It can be written as t ss = l/v s , 
where I and v s are respectively the mean free path and average velocity of the ES during the relaxation process. 
The mean free path will be in turn given by I ~ l/a ss n s , where n s is a density of ESs and a ss is the scattering 
cross-section. 

If we assume the lifetime to be much larger than the inverse frequency of the soliton (t ss ^> 1/lo s ), we can 
resort to the Born approximation to estimate o- ss . As time goes by and the soliton gas rarities due to the process 
of merging, the approximation becomes more and more justified . Within this perturbative treatment and in the 
low velocity limit, we thus find that a ss is simply proportional to the square modulus of J Ui n tdx, where Vint is 
the interaction potential. According to (113) , the latter is proportional to the product of the envelope functions 
of two small-amplitude ESs: Ui„t oc /1/2. Assuming for simplicity that both of them have approximatively 
the same amplitude (f> m ax, it is easy to ascertain from (0) that the overlap integral scales as (f> m ax, similar to 
the energy of a single small-amplitude ES (|l3|). Therefore, the cross-section a ss is proportional to (\> 2 max and, 
according to Eqs. ( p2| ) and (J13|), this implies that a ss ex e. Finally, as n s ex ^/e, see Eq. (pi|), we obtain that 
Zcxe- 15 . 

To estimate v s , we compute first the energy of the moving soliton of Eqs. (pq). In the small-amplitude limit, 
when one can neglect the higher-order and nonlinear dispersive terms, Eqs. ( pq ) admit the exact solution |l|Jl^| 

COSh[(lE - Vgt)lp ma xVO\ 



where V g is a group velocity of the soliton, k = —V g uj max is a soliton quasimomentum and fl — {itpmax 



k 2 /2)/uj max . Using Eqs. (^7j), (|38|), and (|39[), we find the soliton number of quanta A/" s , momentum V s , and 
energy £ s (for uj max = 2): 

A/". = Vwyf , V 8 = -kK=2V g M 8 , £ s = |^ s 3 - 2M s V g 2 = |^ s 3 - ||-. (40) 

These relations show that small-amplitude moving soliton can be indeed considered as a bound state of Af s 
negative-effective-mass quasiparticles (harmonic modes) with a correspondingly positive binding energy equal 

to |a/; 3 (cf. |§). 

If we compare the definition of Af s in (fhj) with expressions (132) and (131) for the energy of the ES (|J) , we 
conclude that Af s oc y/e. Let us now assume that a kind of virial theorem holds for the set of solitons moving 
along the chain with periodic (or fixed) boundary conditions. This amounts to say that (time) averaged kinetic 
and potential energies of such finite dynamical system should be proportional (see, e.g., |lq|). It then follows 
from the definition (WG) of £ s that (V g ) ~ (A/" 2 ) oc e. The validity of such a relation can be physically justified 
by observing that, as a result of the exchange of the number of quanta during the scattering processes, the 
solitons receive the recoil quasimomentum k from free waves which in turn changes the soliton velocities thus 

providing some form of "thermalization" . From all the above follows that v s ~ \ Wo) °^ v 7 ^- 

Finally, we obtain the estimate t ss — l/v s oc e~ 2 and correspondingly t ss 3> l/u) s ~ 1/2, that is fully 
consistent with the initial assumptions. As already mentioned, this scaling is in agreement with the numerical 
observation ||,|(j. 

VII. THE BREATHER 

In the previous Section we explained why the scattering favors the creation of a single localized object. We 
now give a brief account of some of its properties. The h n patterns reported in the left panels of Figs. H and 
illustrate the remarkable fact that a macroscopic fraction of the total energy accumulates on a few sites. For 
instance for the case e = 0.033 and N — 1024, only three sites bear an energy 10.9 corresponding to about 30% 
of the total. As expected from the previous Section, the remaining is distributed along the chain in the form of 
small amplitude (harmonic) oscillations (see again Fig. Eh 

If we observe the breather on a time scale during which its amplitude does not change considerably, we see 
that it randomly alternates between a standing and moving state. Fig. clearly shows how each motion occurs 
with an almost constant velocity and is accompanied by changes of the shape. For a given amplitude, the 
modulus of such a velocity is always the same. We checked that its value is consistent with previous studies 
on exact breather solution: for example from Fig. H we found indeed the velocity to be 0.12(8) in excellent 
agreement with the value 0.120 obtained interpolating the data of Ref. [ pl[ . 

The main issue concerning the velocity is the one of its selection. Although we will not enter into details here, 
we just anticipate that some analytic clue can be achieved by studying the moving solutions of Eq. (0) p3| . In 
particular, we expect the velocity to increase with the amplitude (almost linearly in the large-amplitude limit 
Q4>1nax ^ 1), a fact that is actually consistent with our as well as with other simulations Ipp^l- 

Let us now discuss the power spectra of those strongly localized states. The right panels of Fig. BJ show how 
they differ depending on the amplitude. At smaller energies (upper panels) they resemble more closely an ES 
and the spectrum would have the exponentially decaying prescribed by (|l0|). Upon increasing the energy, Sk 
acquires a sinusoidal component like (|lq ) characteristic of a breather solution. This is simply the consequence 
of the fact that in the continuum approach the two are the only limiting cases of the general solution (M) with 
the amplitude 4> m ax being determined only by the total energy of the lattice. 

Obviously, the existence of the vibrational background affects the structure of the spectra. As the background 
contains a non- negligible energy (the remaining 70% in the case mentioned above), a complete description of 
the spectra requires some information on its features. One consequence that we can draw from Fig. ^ is that 
the naive picture of a single breather moving through a gas of small waves is oversimplified. In fact this would 
imply the spectrum to be a sum of a part like ( |14|) plus an almost flat part. The form of Sk seems rather to 
indicate that only the high-frequency domain of the harmonic spectrum is populated. 



To give a complete description of the relaxation process, we should now consider the final stage of the breather 
destruction (see also Ref. |6fl for a discussion) . This requires understanding how energy is transferred from large 
to small wavenumbers. This issue involves in turn the interaction of a large-amplitude ES with small-amplitude 
plane waves. The latter problem has been studied in Ref. |23|. Here we limit ourselves to observe that it is still 
not clear how to explain why the tipical time scale for such a process should scale again as e~ 2 . An analytic 
argument for that can be actually given only in the very last stage, when the system is sufficiently close to 



equilibrium and relaxation of Fourier modes can be described by usual linear-response theory 24 



VIII. FINITE-TEMPERATURE EFFECTS 

What has been sofar discussed refers to a situation where all the lattice energy is initially fed in the boundary 
mode. This is clearly not realistic for actual experimental conditions, when the system is in contact with a 
thermal bath which populates all the other modes even at (classically) low temperatures. One may therefore 
wonder to what extent the mechanisms described above are effective in localizing the energy. 

In the present section we try to approach the problem by first considering a more general class of initial 
conditions defined by replacing second expression in ( p5| ) by 

u„(0) = {-iTV + a'in , (41) 

where £; is a gaussian random number such that (£;) = and (£z£z+ m ) = <5 m o- In order to assign the initial 
energy per particle to be equal to e, we choose V = 2e/v2e + a 2 and a = er^/2e/(2e + a 2 ). The parameter a 
controls thus the "level of noise" in the initial conditions. Indeed, the limit a 2 <e (exactly the one considered 
in the previous simulations) amounts to feeding all the available energy in the Brillouin-zone-boundary mode 
whereas upon increasing a, we can initially excite also the other Fourier modes of the chain. With our definition, 
equilibrium (i.e. equipartition of energy) is attained for a — > oo (actually when a 2 — > cr^ nax = 2eN 3> e). 

The results presented before can be first extended to the finite a- values by the following reasoning. We showed 
that the initial density of ESs is basically determined by the vibrational energy available in the Brillouin-zone- 
boundary mode. Therefore, one can straightforwardly replace v2e by V in Eq. (|3l|) obtaining 

/ s 1 /3 2e 

7T V 4 V2e + a 1 

As seen in Fig. g, the numerical data are in good agreement with this formula, at least for moderate amplitudes 
a 2 < e/2. This indicates that the mechanism of vibrational energy localization, induced by the modulational 
instability is relatively robust, at least with respect to those type of perturbation in the initial conditions. 

Remarkably, as soon as a overcomes the value e/2, n s seems to be slightly less (around 10 %) than what 
is expected from the purely energetic argument leading to flJ2[ ). Although in this intermediate regime the 
numerical measurements are definitely less accurate due to the large vibrational background, we can at least 
give a qualitative justification to this. In fact, in this case we expect some effective dissipative mechanisms to 
enter into play and some of the localized excitation will be effectively damped out by the interaction with the 
"thermal bath" of the other vibrational modes |24J . 

In order to better clarify this, we performed also some measurements of the power spectrum ( pTJ ) for the 
relatively large values of a (see Fig. R). The results confirm that even for a 2 w e/2 (lower panels) some kind of 
"localization" is still present, although it is hardly detectable by direct measurements. 

To summarize, the creation of localized excitations, induced by the modulational instability, is more and 
more inhibited the closer one starts to the equipartition (thermal equilibrium) state. This is caused both by 
the less amount of energy available in the corresponding short-wavelength mode (as V — > for a 2 ^> e, see 
Eq. (f4l|)) and by the effective dissipation caused by the interaction of small- amplitude ESs with the vibrational 
background. At least in principle, the envelope-function approach could be a useful starting point for a more 
detailed analysis once one would be able to take into account the fluctuating background by casting it in suitable 
stochastic terms. 
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FIG. 1. Modulational instability of the boundary mode for the quartic FPU model with e = 0.0330, N = 1024 The 
left panels are snapshots of the energy density, the right ones the spatial Fourier spectrum of the velocities. The vertical 
lines correspond to the theoretical value of k^ given by Eq. (tCj). 
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FIG. 2. The soliton gas generated after the instability for N = 1024 for e = 0.0033, t = 8000 (upper panels) and 
e — 0.0330, t — 800 (lower panels). The corresponding power spectra Sk of the particles' velocities (right panels) averaged 
over a time interval ~ 100. The dashed lines are the best fit with the formula (p4|). 



0.04 



0.03 



0.02 



0.01 



0.00 











/' 


n ^ 








/ 




n s /E s 






/ 


D- 


0.2 


/Aa 




- 


/ 
/ 






A 


A 


- 0.1 








/ 




0.0 






E 


/t] 










0.00 




0.01 


/ 

/ 
/ 








/ 
/ 


s 






/ 








ur 










/ 








/ 








- 


p.-* 










/ 










/ 










/ 










/ 










< 






I , I 





0.00 



0.01 



0.02 

(3e/2) 1/2 /7i 



0.03 



0.04 



FIG. 3. The maximal initial density n a of ESs as a function of the energy density for N — 1024- . The solid line 
corresponds to formula (|3l|). The inset shows the ratio between the number and energy of the solitons: the horizontal 
line corresponds to Eq. (K3J) with q — 1. 
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FIG. 4. Final breather for different energies e = 0.0033, TV = 256, t = 1.2 10 6 (upper panels) and e = 0.0330, 
TV = 1024, t — 8 10 4 (lower panels). The last spectrum has been averaged over a time 10 4 . Notice the different vertical 
scales. Dashed lines correspond to Eqs. ([lot), (til) and (|lq) with the numerically measured value of 4> m ax and kb — n/3. 
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FIG. 5. The position of the maxima and three successive snapshots of the final breather as a function of time for 
= 0.0330 and TV = 1024. During this observation time the value of the breather energy is practically constant. 
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FIG. 6. The maximal initial density of ESs ns as a function of the parameter a for several values of N ranging 
between 920 and 1024. Solid lines correspond to formula (U2). 



3.00 
2.00 
1.00 



h 



0=0.07 



^\±..^M^u L )M\L,AA,MMU^nL±,Mj',L^,iK 



0.06 



0.03 




0.00 
1000 1 

0.06 



0.03 





0.00 
1000 1 2 fc 3 

FIG. 7. The distribution of energy densities and power spectra of the velocities for e = 0.0330, N = 1024, t = 8 10 
and two different amplitudes of initial perturbation a = 0.07 (upper panels) and 0.12(lower panels). 
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